# 加载相关包
library(openxlsx)
library(geepack)
library(gtsummary)
library(officer)
library(flextable)
# 加载自定义函数
source("./func/gee_run.r")

# 用一个元素为 gee_fit 的二维 list() 绘制所需的表格
summary_out = function (fit_table) {
  # summary_out
  ret = list()
  for (i in 1:length(fit_table[[1]])) {
    tb1 = list()
    for (j in 1:length(fit_table)) {
      tb1[[j]] = tbl_regression(fit_table[[j]][[i]], exponentiate = T, include = c(1))
      print(i)
    }
    ret[[i]] = tbl_merge(tb1, tab_spanner = varsy)
    
  }
  return (ret)
}

###### 读取数据，数据预处理
data_test = read.xlsx("./input/data918.xlsx", 3)
data_anal = data_test[, c(1, 9:13, 14:25, 27, 29:39)]
varsx = colnames(data_test)[c(14:25, 27)]  # 自变量
varsy = colnames(data_test)[c(9:13)]       # 因变量
varsc = colnames(data_test)[c(29:39)]      # 协变量

# 确保变量类型一致：因变量为数值，自变量、协变量为因子
data_anal$id = as.character(data_anal$id)
data_anal[varsy] = lapply(data_anal[varsy], as.numeric)
data_anal[varsx] = lapply(data_anal[varsx], as.factor)
data_anal[varsc] = lapply(data_anal[varsc], as.factor)
# str(data_anal)
######


# 1. 多因素删除缺失值
## 1.1 先删除包含缺失值的行 
tmp_data = data_anal[complete.cases(data_anal), ]
tmp_data[varsx] = lapply(tmp_data[varsx], as.numeric)  # 删除缺失值后，某些因子项的 level 计数可能变为 0，因此需要剔除该 level
tmp_data[varsc] = lapply(tmp_data[varsc], as.numeric)  # 同上
tmp_data[varsx] = lapply(tmp_data[varsx], as.factor)
tmp_data[varsc] = lapply(tmp_data[varsc], as.factor)
## 1.2
fit_table_1 = gee_run(data = tmp_data,
        varsx = varsx,
        varsy = varsy,
        varsc = varsc[-5],
        file = "./output/0919_delete_10cv_del_birthweight_c.xlsx")

s1 = summary_out(fit_table_1)
sect_properties = prop_section(
  page_size = page_size(orient = "landscape",
                        width = 8.3, height = 11.7),
  type = "continuous",
  page_margins = page_mar()
)
save_as_docx(
  values = lapply(s1, as_flex_table),
  path = "./output/921_summary_delete.docx", pr_section = sect_properties
)


## 2.mice
mice_data = read_excel("./output/0919_mice_data.xlsx", 4)
mice_data = mice_data[complete.cases(mice_data[, varsx]), ]
mice_data[varsy] = lapply(mice_data[varsy], as.numeric)
mice_data[varsx] = lapply(mice_data[varsx], as.factor)
mice_data[varsc] = lapply(mice_data[varsc], as.factor)
fit_table_2 = gee_run(mice_data,
        varsx = varsx,
        varsy = varsy,
        varsc = varsc[-5],
        file = "./output/0919_mice_10cv_del_birthweight_c.xlsx")
s2 = summary_out(fit_table_2)
sect_properties = prop_section(
  type = "nextPage",
  page_size = page_size(),
  page_margins = page_mar()
)
save_as_docx(
  values = lapply(s2, as_flex_table),
  path = "./output/921_summary_mice.docx", pr_section = sect_properties
)


# 3.summary_score_qualified
## 3.1 delete
tmp_data = data_anal[, c("id", "score_qualified", varsc)]
tmp_data = unique(tmp_data)
# tmp_data = tmp_data[complete.cases(tmp_data), ]
s3 = tbl_summary(tmp_data[, -1], by = "score_qualified") %>% add_p()

## 3.2 mice
mice_data = read_excel("./output/0919_mice_data.xlsx", 2)[, c("id", "score_qualified", varsc)]
mice_data = unique(mice_data)
s4 = tbl_summary(mice_data[, c("score_qualified", varsc)], by = "score_qualified") %>% add_p()

## 3.3 output
save_as_docx(
  values = lapply(list(s3, s4), as_flex_table),
  path = "./output/921_summary.docx", pr_section = sect_properties
)
